# 	Filename: logdeathsLD_source.R
# 	Description: Source script to produce results based on regressions where "logdeaths" is the dependent variable using the original data with listwise deletion. Called by logdeaths_driver.R.
# 	Author: Barry Hashimoto
# 	Date: August 2019
# 	For: Barry Hashimoto, "Autocratic Consent to International Law: the Case of the International Criminal Court's Jurisdiction," International Organization.
  
############################################################################################################
# Prepare data and match

 data <- NULL
 data <- all
 data$ruleoflaw <- data$vdem.ruleoflaw
 library(caret)
 legaldummies <- as.data.frame(predict(dummyVars(~ as.factor(legaltrad2), data = data), 
 newdata = data))
 names(legaldummies) <- c("civil", "common", "mixed", "islam")
 data <- data.frame(data, legaldummies)
 data$ccode <- factor(data$ccode)
 data$postrome <- ifelse(data$quarter>=154,1,0)  	 
# Identify the OECD DAC States so they can be dropped
  dacmemberlist <- c("Australia", "Austria", "Belgium", "Canada", "Denmark", "Finland", "France", "Germany", 
    "Greece", "Ireland", "Italy", "Japan", "South Korea", "Luxembourg", "Netherlands", "New Zealand", "Norway", 
    "Portugal", "Spain", "Sweden", "Switzerland", "United Kingdom", "United States of America")
 data <- data[ifelse(data$country %in% dacmemberlist, 1, 0)==0,]
 data <- data[data$quarter>=154, ]
 dim(data)

 
############################################################################################################
# Analyze the regime type selected in logdeaths_driver.R

# A function to stop excess verbosity when running matchit().
 quiet <- function(x) { 
   sink(tempfile()) 
   on.exit(sink()) 
   invisible(force(x)) 
  }


# Created matched data sets.

set <- subset(data, subset = democracy == RegimeOption) # Selection of autocracies or democracies made in the driver script.

# Prepare data and match
data.list <- list()
data.list <- with(set, data.frame(leaderexitsoffice, logdeaths, ptssum.mean, RomeStatuteRatification, maleunemp, loggrowth, logpop, logrefugees, rug, eth, rel, logoilrent, logexpmil, logsoldierspercent, capital.scaled, gdp.scaled, ruleoflaw, democracy, common, mixed, islam, pre98conflict, medianIMR, ccode, lcode, year, quarter, quartersinoffice, deaths.trend.single, deaths.trend.country))

# Complete cases.
  data.list <- data.list[complete.cases(data.list), ]

breaks <- 1/2
# Define cutpoints for coarsened exact matching on continuous variables
rug.cuts <- quantile(data.list$rug, seq(0,1,breaks))
logexpmil.cuts <- quantile(data.list$logexpmil, seq(0,1,breaks))
logsoldierspercent.cuts <- quantile(data.list$logsoldierspercent, seq(0,1,breaks))
logpop.cuts <- quantile(data.list$logpop, seq(0,1,breaks))
eth.cuts <- quantile(data.list$eth, seq(0,1,breaks))
rel.cuts <- quantile(data.list$rel, seq(0,1,breaks))
logoilrent.cuts <- quantile(data.list$logoilrent, seq(0,1,breaks))
maleunemp.cuts <- quantile(data.list$maleunemp, seq(0,1,breaks))
loggrowth.cuts <- quantile(data.list$loggrowth, seq(0,1,breaks))
gdp.cuts <- quantile(data.list$gdp.scaled, seq(0,1,breaks))
capital.cuts <- quantile(data.list$capital.scaled, seq(0,1,breaks))
rol.cuts <- quantile(data.list$ruleoflaw, seq(0,1,breaks))
pts.cuts<- quantile(data.list$ptssum.mean, seq(0,1,breaks))
imr.cuts <- quantile(data.list$medianIMR, seq(0,1,breaks))
refugees.cuts <- quantile(data.list$logrefugees, seq(0,1,breaks))

cut.list <- list(rug = rug.cuts, logexpmil = logexpmil.cuts, logsoldierspercent = logsoldierspercent.cuts, logpop = logpop.cuts, eth = eth.cuts, rel = rel.cuts, logoilrent = logoilrent.cuts, maleunemp = maleunemp.cuts, loggrowth = loggrowth.cuts, gdp.scaled = gdp.cuts, capital.scaled = capital.cuts, ruleoflaw = rol.cuts, ptssum.mean = pts.cuts, medianIMR = imr.cuts, logrefugees = refugees.cuts)

match.list <- quiet(matchit(RomeStatuteRatification ~ logexpmil + logsoldierspercent + logpop + eth + rel + rug + logrefugees + logoilrent + loggrowth + maleunemp + gdp.scaled + ruleoflaw + capital.scaled + common + mixed + islam + pre98conflict + medianIMR + ptssum.mean, data = data.list, method = "cem" , cutpoints = cut.list))

#summary(match.list)

# No controls, shared trend
shared.nocontrol <- lm(formula = logdeaths ~ RomeStatuteRatification + deaths.trend.single + as.factor(ccode) + as.factor(year), data = set)

# Controls, shared trend
shared.control <- lm(formula = logdeaths ~ RomeStatuteRatification + logexpmil + logsoldierspercent + logpop + eth + rel + rug + logrefugees + logoilrent + loggrowth + maleunemp + gdp.scaled + ruleoflaw + capital.scaled +  deaths.trend.single + as.factor(ccode) + as.factor(year), data = set)

# No controls, unit-specific trend
country.nocontrol <- lm(formula = logdeaths ~ RomeStatuteRatification + deaths.trend.country + as.factor(ccode) + as.factor(year), data = set)

# Controls, unit-specific trends. 
country.control <- lm(formula = logdeaths ~ RomeStatuteRatification + logexpmil + logsoldierspercent + logpop + eth + rel + rug + logrefugees + logoilrent + loggrowth + maleunemp + gdp.scaled + ruleoflaw + capital.scaled +  deaths.trend.country + as.factor(ccode) + as.factor(year), data = set)

# Matched data, Battle deaths, shared trends
matched.shared <- lm(formula = logdeaths ~ RomeStatuteRatification + logexpmil + logsoldierspercent + logpop + eth + rel + rug + logrefugees + logoilrent + loggrowth + maleunemp + gdp.scaled + ruleoflaw + capital.scaled + common + mixed + islam + pre98conflict + medianIMR + ptssum.mean + deaths.trend.single, data = match.data(match.list))

# Matched data, Battle deaths, unit-specific trends
matched.country <- lm(formula = logdeaths ~ RomeStatuteRatification + logexpmil + logsoldierspercent + logpop + eth + rel + rug + logrefugees + logoilrent + loggrowth + maleunemp + gdp.scaled + ruleoflaw + capital.scaled + common + mixed + islam + pre98conflict + medianIMR + ptssum.mean + deaths.trend.country, data = match.data(match.list))

satt.battles <- zelig(formula = logdeaths ~ RomeStatuteRatification + logexpmil + logsoldierspercent + logpop + eth + rel + rug + logrefugees + logoilrent + loggrowth + maleunemp + gdp.scaled + ruleoflaw + capital.scaled + common + mixed + islam + pre98conflict + medianIMR + ptssum.mean + deaths.trend.country, model = "ls", data = match.data(match.list), cite = F) %>% ATT(treatment = "RomeStatuteRatification", treated = 1, num = 20000) %>% get_qi(qi="ATT", xvalue = "TE")


  writeLines("")
  writeLines("# Print 95% interval plus 25th and 7th percentiles (top) and median absolute deviation (bottom) of the SATT of RomeStatuteRatification on logdeaths for states of the selected regime type. Unimputed data with listwise deletion.")
  print(round(quantile(satt.battles, c(0.025, .25,.5,.75,.975)), 3))
  print(round(mad(satt.battles), 3))


# A function to get the adjusted R squared from lm()
 get.r.lm <- function(x){return(summary(x)$adj.r.squared)}
  
 output.deaths <- function(x){
	mat <- rbind(summary(x)$coefficients[2,c(1:2, 4)], 
	             c(get.r.lm(x), NA, NA),
	             c(length(x$residuals), NA, NA)
	             )
	mat[c(1:2),] <- signif(mat[c(1:2),], 3)            
	rownames(mat) <- c("Coef. of ICC Jurisdiction", "Adj. R-squared", "N")
	return(mat)
	}

 col1 <- rbind(output.deaths(shared.nocontrol), output.deaths(shared.control), output.deaths(matched.shared))
 col2 <- rbind(output.deaths(country.nocontrol), output.deaths(country.control), output.deaths(matched.country))

 deathstable <- cbind(col1, col2)
 colnames(deathstable) <- rep(c("Est.", "SE", "p"), 2)

  writeLines("")
  writeLines("# Print models of logdeaths for the selected regime type using the unimputed data with listwise deletion. First column has models with shard trends, while second column has models with country-specific trends. First row prints results for models with country and year fixed effects but WITHOUT control variables. Second row prints results for models with the same fixed effects and WITH control variables. Third row prints results for models using the matched data sets WITH control variables.")

  print(round(deathstable, 3))

# library(memisc)
# toLatex(deathstable)

  writeLines("")
  writeLines("Total number of rows in the set of observations for this regime type since 1998-Q3 after listwise deletion.")
  print(dim(data.list)[1])
 
  writeLines("")
  writeLines("Total number of states of this regime type since 1998-Q3 after listwise deletion.")
  print(length(unique(data.list$ccode)))

  writeLines("")
  writeLines("Total number of unique leader IDs in states of this regime type since 1998-Q3 after listwise deletion.")
  print(length(unique(data.list$lcode)))
